{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# DP mixture of Gaussians as used in Anglican\n",
    "using Turing, Distributions\n",
    "data = [1.0, 1.1, 1.2, -10, -15, -20, 0.01, 0.1, 0.05, 0]\n",
    "N    = length(data)\n",
    "\n",
    "@model dpmixa begin\n",
    "  @assume ϕ ~ IArray(Normal(0, 10)) # Infinite Array\n",
    "\n",
    "  urn = PolyaUrn(1.72)\n",
    "    μ, classes = tzeros(N), tzeros(Int, N)\n",
    "  for i in 1:N\n",
    "    classes[i]  = randclass(urn)\n",
    "    μ[i]  = ϕ[classes[i]]\n",
    "    @observe data[i] ~ Normal(μ[i], 1)\n",
    "  end\n",
    "  K = length(unique(classes))\n",
    "\n",
    "  @predict K μ\n",
    "end\n",
    "\n",
    "# Collect 50 samples from Particle Gibbs sampler\n",
    "@time resulta  = sample(dpmixa, PG(20, 50))\n",
    "\n",
    "println(resulta[:K])\n",
    "#println(resulta[:μ])\n",
    "\n",
    "println(\"Model evidence (: $(resulta[:logevidence])\")\n",
    "println(\"$(mean(resulta[:K])), $(mean(resulta[:μ]))\")\n",
    "\n",
    "\n",
    "@time resulta2  = sample(dpmixa, SMC(10000))\n",
    "println(\"$(mean(resultb2[:K])), $(mean(resultb2[:μ]))\")\n",
    "println(\"Model evidence (log): $(resulta2[:logevidence])\")\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.5",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
